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Human perception, cognition, and action are supported by a complex network 
of interconnected brain regions. There is an increasing interest in measuring and 
characterizing these networks as a function of time and frequency, and inter-areal phase 
locking is often used to reveal these networks. This measure assesses the consistency 
of phase angles between the electrophysiological activity in two areas at a specific 
time and frequency. Non-invasively, the signals from which phase locking is computed 
can be measured with magnetoencephalography (MEG) and electroencephalography 
(EEG). However, due to the lack of spatial specificity of reconstructed source signals in 
MEG and EEG, inter-areal phase locking may be confounded by false positives resulting 
from crosstalk. Traditional phase locking estimates assume that no phase locking exists 
when the distribution of phase angles is uniform. However, this conjecture is not true 
when crosstalk is present. We propose a novel method to improve the reliability of 
the phase-locking measure by sampling phase angles from a baseline, such as from a 
prestimulus period or from resting-state data, and by contrasting this distribution against 
one observed during the time period of interest. 
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INTRODUCTION 

Indices of consistent phase differences at a particular time and 
frequency are commonly used to assess coherence relationships 
between two signals. It is calculated by first computing the 
time-frequency representations of a pair of time-varying sig- 
nals. This may be accomplished through a variety of means, 
including wavelet filtering and by applying the Hilbert trans- 
form. Thereafter, the phase difference at the corresponding time 
and frequency points between the two time courses is computed 
on a trial-by-trial basis to test for repeatability. In magnetoen- 
cephalography (MEG) and electroencephalography (EEG), phase 
locking is evaluated either between two sensor signals, between 
two estimated source time courses, or between one sensor or 
source time course and an external reference signal such as the 
electromyogram (EMG) (Tallon-Baudry et al., 2001; Lin et al., 
2004; Schoffelen and Gross, 2009). 

Constant phase differences between activation time courses 
can be attributed to functional connectivity. In particular, two 
areas whose estimated source time courses have a consistent phase 
difference are likely communicate with each other through some 
pathway or may be jointly affected by a third source. 

If the phase relationship is not fixed, the difference of the 
phase angles will be random with no bias toward a specific angle. 
Therefore, under the null hypothesis of no phase consistency, the 
angles will have a uniform circular distribution. However, in real- 
ity, various unavoidable confounds in data processing may lead to 
false consistent phase differences. 



Non-invasive recordings of human electrophysiological activ- 
ity can be obtained with EEG and MEG. However, mapping the 
MEG/EEG sensor signals back onto the cortex is an ill-posed 
inverse problem (Hamalainen et al., 1993) and, therefore, appro- 
priate constraints and regularization need to be applied to render 
the solution unique and stable. One common solution known 
as the minimum-norm estimate (MNE), finds a smooth, dis- 
tributed current distribution with minimum L2 norm constraint 
to reproduce the measured data (Hamalainen and Ilmoniemi, 
1994). With help of high-resolution MRI data, it is possible to 
reconstruct the geometry of the individual cortical surface (Dale 
et al, 1999; Fischl et al, 1999a,b, 2001) and use this information 
to create location and orientation constraints for source modeling 
(Dale and Sereno, 1993; Dale et al., 2000; Lin et al, 2006). 

However, due to its distributed nature, the MNE will exhibit 
strong false spatial correlations due to spread and cross talk in 
the estimates, which often lead to such false positives (Schoffelen 
and Gross, 2009). Beamforming methods (Van Veen et al., 1997; 
Robinson and Vrba, 1999; Hui et al., 2010; Rana et al., 2011), 
designed for spatial filtering of data, help to reduce correlations, 
but will still tend to map the same sensors, with variable weights, 
onto multiple ROIs, and this can also lead to false positives in 
phase locking. 

With direct measurements of neuronal activity through single- 
unit, multi-unit, or local field potential (LFP) recordings spaced 
apart to avoid strong correlations, phase locking methods can 
be applied without concerns of confounds present in MEG/EEG 
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source modeling. However, the resulting data will only allow for 
inferring phase locking amongst a small set of neurons instead 
of between functional brain areas. Thus, invasive recordings only 
solve the issue of crosstalk but cannot be efficiently performed 
across the cortex. In addition, direct invasive recordings are very 
rarely conducted in humans and, when they are, are used in clin- 
ical cases such as in epileptic patients (Gloor, 1975; Keene et al, 
2000). 

In this methodological paper, we propose a novel, fast method 
with which the distribution of phase angle differences between 
two areas are compared to an arbitrary null distribution, e.g., to 
the distribution of phase angles in the absence of the stimulus, via 
a non-parametric phase angle method with asymptotic bounds, 
producing significance levels with minimal computational cost. 
We will describe the traditional phase locking method and com- 
pare it to our new method, Uniform Scores Test for Phase Locking 
(USTPL). In addition, we shall discuss two recently introduced 
methods that also provide a means for baseline correction: the 
Phase Lag Index (PLI) and the Phase Bifurcation Index (PBI). 
We will show how non-parametric statistical methods applied 
directly to the phase angle differences result in an accurate sig- 
nificance test as opposed to using a parametric test that falsely 
assumes a uniform null distribution. 

METHODS 

The first section here discusses the nature of the simulated data 
and how data is mapped onto the brain. The second section elab- 
orates on the four methods that are the focus of the comparison: 
traditional phase locking, our proposed method USTPL, as well 
as PLI and PBI. 

DATA SIMULATION 

Brain surface and region of interest (ROIj selection for data 
simulation 

The brain surface used to produce simulations is from a healthy 
male subject, age 23 at time of collection. Structural MRI scans 
were acquired using an 8-channel phase-array head coil in a 3T 



scanner (Siemens-Trio, Erlagen, Germany). Freesurfer software 
(http://surfer.nmr.mgh.harvard.edu) (Fischl et al, 2002, 2004) 
was used for cortical surface reconstruction and parcellation. 
ROIs (Figure 1) were chosen based on clusters of significant acti- 
vation during a visual search task with a facilitatory auditory cue 
adapted from (Vaina et al, 2010; Calabro et al, 201 1). 

Data simulation 

Simulated data consisted of 20, 50, or 80 trials with phase locking 
between all ROI pairs and 20, 50, or 80 trials with no phase lock- 
ing. The trials with no phase locking serve as a means for baseline 
comparison for the USTPL [see section "Uniform Scores Test for 
Phase Locking (USTPL)] and PBI (see section "Phase Bifurcation 
Index") methods. Phase locking was simulated through the fol- 
lowing steps. First, we randomly selected, from a uniform circular 
distribution, a fixed phase angle for each ROI to ensure that all 
ROIs will be phase locked. Next, in order to perturb the phase 
angle slightly between trials, we added a zero-centered von-Mises- 
distributed random angle (Mardia and Jupp, 2000) with a varying 
concentration parameter k = (20, 50). We define the von-Mises 
distribution as follows: 



me, k) 



cos(x— 8) 



2jtJo0c) 



(1) 



where 6 is the phase angle at which the distribution is centered, 
k is a measure of concentration of the distribution around 0, and 
7o is the modified Bessel function of order 0. 

We then produced a 40-Hz sinusoidal time course in each ROI 
with a phase shift equal to the perturbed phase angle associated 
with that ROI. This was repeated for each trial and we thus obtain 
simulated source-space waveforms xjt(f) for each trial k. 

To accurately model the MEG data measurements as well as 
to induce cross-talk, we performed a multi-step process summa- 
rized in Figure 2. First, simulated data in each ROI, generated in 
the cortical source space was mapped to the sensor space using the 
forward operator (gain matrix) G, which was computed using a 
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FIGURE 1 | Regions of interest chosen from which to generate simulation data. ROIs are adapted from Vaina et al. (2010); Calabro et al. (2011). There is no 
overlap between ROIs to avoid any issues in computing phase locking connectivity between them. 
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FIGURE 2 | Flowchart of data simulation. A dotted line in the middle 
separates source and sensor space computations. A random perturbation 
phase angle is added to the fixed phase angle at each ROI to induce a 
phase shift on the sinusoidal time course of that ROI. The time course is 



mapped to the sensor space, sensor noise is added and the simulated 
noisy signals are mapped back to the source space using the MNE linear 
inverse operator. In trials without phase locking, the fixed ROI phase angle 
was also randomized. 



single-compartment boundary-element model with the shape of 
the inner skull surface extracted from the MRI data of the sub- 
ject (Hamalainen and Sarvas, 1989). Next, we added in sensor 
noise n^t) to each trial. This was obtained from the pre-stimulus 
baseline data of the experiment described in Vaina et al. (2010). 
We first computed an estimate for the noise covariance matrix 
from these data and then used this matrix to obtain spatially 
colored noise from Gaussian noise with unit variance and zero 
mean, independent across sensors. As a result we obtain the noisy 
sensor-space signals: 

y k (t) = Gx k (t) + n k (t) (2) 

As an inverse solution we employed the cortically constrained, 
depth-weighted L2 MNE (Hamalainen and Ilmoniemi, 1994; 
Dale et al., 2000; Lin et al., 2006). These computations were per- 
formed using the MNE software (http://www.nmr.mgh. harvard. 
edu/martinos/userlnfo/data/sofMNE.php), which amounted to 
multiplying yt(f ) with a linear inverse operator W: 

x' k (t) = Wy k (t) = W(Gx fc (0 + n k (t)) (3) 

In the computation of W, we constrained the source orientations 
to be normal to the cortex and used the same noise-covariance 
matrix as for generating a.k(t). For each we then computed an 
average MNE waveform for each ROI. To avoid signal cancella- 
tion in the averaging, we calculated the principal direction of the 
cortical surface normals within each ROI and inverted the sign of 
the waveform at the given vertex if the surface normal at this ver- 
tex pointed to a direction opposite to the principal surface normal 
direction. 

In trials with no phase locking, the phase angle associated with 
each ROI was randomized on a trial-by-trial basis. This ensured 
an equal signal amplitude without any phase-synchrony. 

Wavelet filtering for phase locking 

We extracted phase time courses in each ROI for a particu- 
lar frequency by wavelet filtering. The trial-by-trial ROI time 



courses were decomposed into complex time-frequency coeffi- 
cients through the Morlet wavelet transform (Lin et al., 2004) 
with envelope bandwidth of 1/5 octaves. The time courses in each 
ROI were filtered with a Morlet wavelet at 40 Hz at 0.5 s into the 
stimulus to obtain a complex time and frequency filter coefficient 
for each trial. 

Computation of receiver operating characteristic curves 

A Receiver Operating Characteristic (ROC) curve is a representa- 
tion of the error of classification over varying thresholds (Green 
and Swets, 1974). It is a plot of the true positive rate (TPR) vs. 
the false positive rate (FPR) of a detector, which in this case is 
the detection of significant phase locking. This method is ideal to 
compare the four phase locking methods since the PLB and PLI 
statistics will not provide a significance level directly. 

For each set of parameters for each method, we repeat the data 
simulation of section "Data Simulation" 20 times to obtain a set of 
statistics from which to construct the ROC curve. We computed 
one set of 10 simulations to produce positives, with a set of phase- 
locked trials. Another set of 10 simulations produced negatives, 
with no phase-locked trials. For each set of parameters and each 
iteration of the methodology from section "Data Simulation", we 
obtained a single statistic at 0.5 s into the sample (see section 
"Wavelet Filtering for Phase Locking") from every ROI compari- 
son. Since there are 22 ROIs, each simulation iteration produced 
231 comparisons. Combining across all simulations, we obtained 
2310 positives and 2310 negatives total. 

We sorted the statistics from all iterations by their statistical 
value and used each as a threshold, such that a value that is equal 
to or more significant than the current threshold was labeled a 
positive and ones that are lower were labeled a negative. Thus, we 
computed two values that characterize the performance of each 
phase locking detection method: The TPR, which is the rate at 
which positives were correctly classified as positives, and the FPR, 
which is the rate at which negatives were incorrectly classified as 
positives. The TPR and FPR were computed as TPR = TP/P and 
FPR = FP/N, respectively. Here, TP is the number of positives 
accurately evaluated as significant and P is the total number of 
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true positives while FP is the number of false positives incorrectly 
classified as significant and N is the total number of true nega- 
tives. For every threshold, we obtained a single point on the ROC 
curve, using the TPR and FPR. Thereafter, we sorted the values 
according to the threshold and connected the ordered points to 
obtain the final ROC curve. 

We characterized each ROC curve with the area under the 
curve (AUC). For a completely random assignment to positive 
and negative detection, the ROC curve is a line TPR = FPR result- 
ing in AUC = 0.5. When TPR > FPR, the AUC increases reaching 
AUC = 1 when TPR = 1 and FPR = 0. 

Error bounds were computed from the above procedure com- 
puting the AUC over 5 sets of 4 iterations (using the 20 trials) to 
generate 5 separate ROC curves for each set of parameters. The 
mean is reported as the AUC value and the standard deviation is 
reported as the interval. 

Computation of error rates 

When a statistical significance is readily available, in the case of 
PLV and USTPL, we may compute an error rate. Using a signif- 
icance level at 0.05 with and without Bonferroni correction, we 
computed the error as: 



Error = (TP + TN)/(P + N) 



(4) 



This is the number of correctly assigned positive and negative 
phase locking detections divided by the total number of tests. 

PHASE LOCKING METHODS 

We illustrate with an example the shortcomings of the traditional 
phase locking method, which tests for similarity of trial-by-trial 
signal phase angle differences between two sources. If the differ- 
ence of the phase angles is uniformly distributed across trials, 
the traditional method assumes that there is no phase locking, 
since there is no coherence between the phase angles. However, 
crosstalk may induce a spurious coherence of phase angles. Our 
method extends the traditional phase locking method by test- 
ing the phase angle differences tested against an empirical null 
distribution that we can sample from the pre-stimulus interval. 

Traditional phase locking 

The concept of phase locking lies in the idea of phase as a shift in 
a signal. Figure 3 illustrates this using three pairs of time courses. 
Note that the signals in ROI 1 and ROI 2 occur at random times 
but their relative phase lag is consistent. We can use circular 
statistics to assess whether this lag, across trials, is consistent. 

To compute the phase lag at different frequencies, we uti- 
lize wavelet filtering that produces phase estimates for each time 
point and for desired frequencies, see section "Wavelet Filtering 
for Phase Locking". We can then calculate the trial-by-trial phase 
differences between the two ROIs as a function of time and 
frequency. 

The traditional phase locking methodology assumes that the 
two time courses, if not phase locked are completely independent, 
and thus the phase angle differences will be uniformly random. 
Therefore, one can apply a test of uniformity to assess whether 
the distribution of phase angle differences is indeed uniform. The 
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FIGURE 3 | Illustrative example of phase-locked signals. The red line is a 
reference relative to the peak of the signals. The green line illustrates a 
consistent lag between ROI 1 and ROI 2 across trials. 



most common uniformity test is the Rayleigh test for circular uni- 
formity (Jervis et al., 1983; Tallon-Baudry et al., 1996; Lin et al, 
2004). The Rayleigh test considers each phase angle as a vec- 
tor with unit length and angle equal to the test angle. First, the 
phase angle difference between two ROIs is computed for each 
trial and the corresponding unit vector is found. Then, the vectors 
are averaged across trials to produce a single averaged vector. The 
magnitude of the average vector is the Phase Locking Value (PLV). 
This value ranges from zero to one, corresponding to non-existent 
and complete phase locking, respectively. When the number of 
phase angles averaged is large [n > 50), the significance of the 
phase locking value can be approximated (Lin et al, 2004) by: 



Pplv ~ exp(-PLV). 



(5) 



Crosstalk and phase locking 

While the traditional phase locking method works well for two 
sources whose time courses can be determined without inter- 
ference, the cortical source estimates in neighboring regions, 
computed from MEG sensor signals, are prone to crosstalk (Liu 
et al, 2002) leading to false-positive phase locking detection as 
discussed in section "Data Simulation." 

To demonstrate this, we drew the phase angle differences as 
samples from a von Mises distribution. We used this distribution 
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since it allowed us to control sampling from a unimodal dis- 
tribution with a specified mean and variance, through 9 and 
k, respectively. As k increases, the probability of sampling 
closer to 9 increases. The von-Mises distribution approximates 
the circular normal distribution, which is analogous to the 
normal distribution for circular statistics (Mardia and Jupp, 
2000). 

We generated two sets of 20 samples from a von Mises distri- 
bution (9 = 0, k = 4) to simulate sampling from a pre-stimulus 
interval and a post-stimulus interval where there is no interaction 
between the two ROIs, see Figure 4. A Rayleigh test applied to the 
PLV calculated from the stimulus samples showed that this PLV 
is significant (a < 0.05). However, the distribution in the pre- 
stimulus interval is similar to the post-stimulus interval, implying 
that there is no change in phase locking between the pre- and 
post-stimulus intervals. Therefore, it is necessary to perform a sta- 
tistical test that compares the two distributions in a fast, efficient 
manner. 

While one could assume that it should be sufficient to compare 
the average vector magnitudes via permutation tests or bootstrap 
resampling, these procedures computationally taxing procedures 
and may also lose useful phase information, as illustrated in 
Figure 5A. The PLV's are similar but one can see that there is a 
phase shift between the two. Specifically, the pre-stimulus phase 
angle differences are distributed near zero whereas the post- 
stimulus phase angle samples are close to jt/4. This suggests the 
possibility of a communication lag, or a potential change in the 
delay between information passed through ROI 1 and ROI 2 from 
the pre- to post-stimulus period. To account for these differences, 



we propose a statistical test that compares two distributions of the 
phase angles. 

Uniform scores test for phase locking (USTPL) 

We propose the Uniform Scores Test (UST) (Mardia and Jupp, 
2000) as a method to compare two phase angle distributions. 
In this test, the samples of phase angles from the two distribu- 
tions are linearly ranked jointly so the rank 1 data point will be 
closest to zero degrees in the positive direction, the rank 2 data 
point will be the second closest, while the last data point will be 
furthest from 0, or closest to 2jc. The ranked phase angles are 
then evenly positioned around the unit circle; the corresponding 
polar angles are then called the circular rank statistics (Mardia 
and Jupp, 2000). This is illustrated in Figure 6. 

Using the ranked phase angles, we next compute a statistic 
similar to the phase locking value: 

%=\Z«»tf>) +(l>ff) (6) 

where is the size of phase angle population considered, k, (in 
this case, we can choose the post-stimulus population) and P-*' 
is the circular rank statistic of the ith element of population k. To 
take advantage of asymptotic bounds for large sample sizes, we 
normalize Rk to create a new statistic: 

K = (7) 
nyni 

As the total number of samples grows large (n\ + «2 = n > 40), 
the null distribution approaches a x 2 distribution with two 
degrees of freedom (Mardia and Jupp, 2000). 

For demonstrative purposes, we apply USTPL's ranking 
method in Figure 5A to obtain Figure 5B. The redistribution of 
angles clearly illustrates the separation of the two distributions 
by phase, and thus the difference between the prestimulus and 
poststimulus intervals. 

USTPL for comparing stimulus conditions 

Since the UST is a comparative statistical test, it can be used to 
compare any pair of phase angle populations, e.g., two stimulus 
conditions A and B. Instead of employing UST between the post- 
stimulus and pre-stimulus periods, we can compute it between 
the conditions A and B. One drawback is that, although this 
test will allow detection of a significant phase locking difference 
between A and B, the test alone cannot tell which is stronger 
or whether the phase difference has changed between conditions 
A and B. 

Phase lag index 

A recently proposed approach to solve the cross-talk problem uses 
a statistic known as the PLI (Stam et al, 2007; Vinck et al., 2011). 
In this method, it is assumed that a common source of noise or 
cross-talk between two sources is associated with a phase angle 
of either 0 or jt, the latter corresponding to a change in polarity 
only. PLI, therefore, sets out to measures the asymmetry between 
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FIGURE 4 | Simulated phase angle differences for a poststimulus 
period (Stim) and a prestimulus period (Prestim), both generated from 
a von Mises distribution with 8 equal to zero and k equal to 4. The 

vectors represent the vector averages of the data points corresponding to 
each population of phase angle differences. 
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FIGURE 5 1 (A) Simulated phase angle differences from von Mises their corresponding circular rank statistic. The vectors represent the 

distributions for a post-stimulus period (8=u/4, k = 4) and a vector averages of the data points corresponding to each population 

pre-stimulus period (8 = 0, k = 4). (B) Same data redistributed by of circular ranks. 
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FIGURE 6 | An example of circular rank statistics. The circle on the 
left shows points before circular ranking whereas on the right the 
points are ranked and repositioned. The numbers indicate the rankings 
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of the points. The repositioned points are equally spaced around the 
unit circle and the order of the classification of the red and green 
points is preserved. 



the number of phase angles whose sines are positive and negative. 
The PLI statistic is defined as: 



PLI = |(sign[sin9])| 



(8) 



where 9 is the phase angle difference. If there is a bias of phase 
angle differences toward either side, the PLI will deviate from 
zero, which indicates significance. In our simulations (see sec- 
tion Crosstalk and Phase Locking"), we did not discriminate 



between leading or lagging phases, therefore, in our analysis, we 
do not differentiate between a resulting mean sign that is above 
or below zero. Therefore, we take the magnitude of the mean sign 
as the PLI. 

Vinck et al. (2011) introduce a weighting term to form a 
weighted PLI (WPLI): 



WPLI : 



|3 {X} | sign [sine ])| 

<PTOI> 



(9) 



Frontiers in Neuroinformatics 



www.frontiersin.org 



February 2013 | Volume 7 | Article 3 | 6 



Rana et al. 



Comparative statistical phase locking test 



with 3 {X} = A sin9, where A is the magnitude of the signal. 
Thus, for signals that are small, the PLI score will not be as 
strongly affected. In our simulation the signal amplitudes are held 
constant between the baseline and the stimulus, and differences 
between PLI and WPLI would occur only due to the presence of 
noise. Since the signal amplitude is held constant in our simu- 
lations and that amplitude measures are not incorporated into 
other methods, we use Stam's formulation of PLI to compare 
against the other methods. 

Phase bifurcation index 

We will also consider another method to that uses a similar con- 
cept of a phase angle distribution comparison, is the PBI proposed 
by Busch et al. (2009) defined as: 

PBI = (PLVi - PLV aU ) x (PLV 2 - PLV aU ) (10) 

where PLVi, PLV2, and PLV a ii are the phase locking values from 
distribution 1, 2, and the combined population, respectively. The 
phase locking value of the individual distributions (PLVi, PLV2) 
are calculated by applying the method in section "Traditional 
Phase Locking" to the subset of trials from which the distribu- 
tion (either 1 or 2, such as the prestimulus or poststimulus) is 
sampled. To compute PLV a u, the method should be applied to the 
grouped trials from 1 and 2. Each of the differences in Equation 
(10) compares each population to the combined one. Due to the 
multiplication of the two differences, the PBI statistic will be sig- 
nificant if both populations are significantly different from the 
combined population. 

Obtaining statistical significance measures 

In the above methodologies, only PLV and USTPL provide a 
direct method for obtaining a statistical significance level. In order 
to obtain statistical significance with PLI or PBI, non-parametric 
methods may be used. Since we will be performing a comparison 
between two populations, a simple permutation test should act 
as an appropriate means to compute a difference. For PLI, we can 
do so by computing the PBI during the stimulus and nonstimulus 
interval over a number of iterations while permuting between the 
two. However, this means that the statistic has to be computed a 
large number of times. A conservative number would be to per- 
form 1000 permutations. For PBI, samples for PLVi and PLV2 can 
be permuted between each other, but this requires recomputing 
PLV for each 1000 times. 

Table 1 describes the time required to compute the signifi- 
cance of a single time -frequency datapoint tested on a Macbook 
Pro (version 8,1, Late 2011) with 100 trials. As we sample a 
larger number of time-frequency points, the computation time 
increases linearly, and thus multiplicatively if computed across 
a full grid of time and frequency coordinates as we increase the 
samples on either axis. It is clear that for such scenarios, it would 
be computationally expensive to compute the bootstrap for PLI 
or PBI. 

RESULTS 

We first discuss the performance of our proposed method, 
USTPL, compared to PLI, PBI, and PLV in detecting phase 



Table 1 | Time for computing a single significance statistic. 



Method Computation time (s) 



USTPL 0.01538± 0.00014 

PLI 0.01520 ±0.00014 

PBI 15.2181 ±0.1416 

PLV 15.2199± 0.1415 



This table illustrates the difference in time for computing the PLI and 
PBI statistics, requiring expensive non-parametric techniques (in this case 
permutation testing). USTPL and PLV can be computed in milliseconds whereas 
it takes seconds for PLI or PBI for obtaining a single significance statistic. 
Calculations were performed over 50 trials for each the stimulus and prestimulus 
periods on a Macbook Pro (version 8, 1, Late 2011). The times are reported as a 
mean ± the standard error. 



locking across trials at a particular time slice with help of ROC 
curves (Green and Swets, 1974). Second, we compare the error 
rates of USTPL and PLV for a given statistical significance level. 
Since PLI and PBI do not produce significance levels directly, we 
do not consider them in our comparison. 

DETECTION PERFORMANCE 

The comparison of the ROCs of the four methods shown in 
Figure 7 indicates that, overall, USTPL (shown in blue) has a 
higher true-positive rate than the other three methods in the con- 
dition of a mean noise level equal to 0.1% of the signal level, 
random angle concentration parameter k = 50, number of trials 
n = 20. Although USTPL's ROC curve approaches a TPR of 50% 
as the FPR approaches 50%, this is not likely to be inherent to the 




False Positive Rate 

FIGURE 7 I ROC curve comparing the four statistics at a mean noise 
level equal to approximately 0.1% of the signal level, k = 50, number 
of trials (n) = 20. USTPL (blue) has more detection power than PLV (cyan). 
PLI (red), and PBI (green). 
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Table 2 | Performance of all three methods for various noise levels, number of trials (n), and concentration factors (ic). 



Noise Iv: n = 50, k = 50 n = 20, k = 50 





0.10% 


1.00% 


0.10% 


1.00% 


USTPL 


0.5917 ±0.0839 


0.4962 ±0.0117 


0.5900 ±0.0387 


0.4974 ±0.0233 


PLI 


0.6460 ± 0.0744 


0.5119 ±0.0203 


0.5025 ±0.0271 


0.4983 ±0.0139 


PB 


0.5371 ±0.0233 


0.4938 ±0.0128 


0.5246 ±0.0282 


0.5144 ±0.0174 


PLV 


0.5032 ±0.0164 


0.4966 ±0.0126 


0.5198± 0.0229 


0.4840 ±0.0037 


Noise Iv: 


n = 80, k = 50 




n = 50, k = 20 






0.10% 


1.00% 


0.10% 


1.00% 


USTPL 


0.5989 ±0.0639 


0.5052 ±0.0067 


0.4996 ±0.01 87 


0.5051 ±0.0156 


PLI 


0.6693 ±0.0663 


0.4985 ±0.0118 


0.4991 ±0.0115 


0.5106± 0.0113 


PB 


0.5190 ±0.0072 


0.5100 ±0.0108 


0.5101 ±0.0125 


0.5048 ±0.0126 


PLV 


0.5038 ±0.0173 


0.4914 ±0.0163 


0.5051 ±0.0199 


0.5116 ±0.0146 


The numbers 


are reported as the mean ± 1 standard deviation 


over 5 simulation repetitions 


(section "Computation of Receiver Operating Characteristic 



Curves "). 
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FIGURE 8 | Comparison of error rates at p = 0.05 with and without 
Bonferroni correction between PLV and USTPL methods. Dotted lines 
indicate error rates at Bonferroni-corrected significance values. USTPL is 
significantly lower than PLV at all tested noise levels. 



method but is unique to this particular simulation at these param- 
eter values. Table 2 summarizes the AUC values across different 
conditions. 

As k decreases, we find a significant drop in performance of all 
methods. As number of trials (n) decreases significance strength 
decreases most for PLI. Signal level has a significant effect on 
detection performance as well. 

SIGNIFICANCE LEVEL ERROR RATE 

Figure 8 shows that for a significance level threshold p = 0.05, 
we obtain lower error rates when using USTPL compared to PLV. 
Interestingly, the Bonferroni-corrected significance levels cause 
an increase in the error rates. This is understandable because the 
correction will increase the significance of the threshold, and thus 
true positives with low-significances will be incorrectly classified 
as false. However, this correction is appropriate in experimental 
situations when true positives are fewer and false positives are 
essential to be removed. 

It is interesting that in this simulation, the traditional phase 
locking method has nearly no detection power. However, this does 
not mean it has no use in MEG. Due to the closeness of ROIs on 
the brain in this simulation, the cross-talk is significant. Thus, 
leakage into neighboring ROIs is likely to induce a detectable 
bias to the phase differences between these neighboring ROIs. In 
addition, our simulations were conducted at a single frequency. 
Therefore, all sources, even in the baseline, will generate signals 
at this frequency and thus induce cross-talk. In a more realistic 
situation, there would be activations at different frequencies and 
thus there will be a smaller set of areas where activation occurs at 
a certain frequency. 

DISCUSSION 

This paper provides a novel statistical method to compare the 
phase angle differences between two ROIs to an empirical null dis- 
tribution using the uniform scores test, producing a statistical sig- 
nificance value. We compared this new method (USTPL) to three 



existing methods, PLV, PLI, and PBI. Unlike PLI or PBI, USTPL 
does not require computationally expensive non-parametric sta- 
tistical methods significance levels. Thus this method can be easily 
used for testing many frequencies and time points between a large 
number of ROIs, or for computing statistics across the cortical 
surface as in Lin et al. (2004) for PLV. 

In our simulation, we were able to detect phase locking when 
the concentration parameter was high and noise was low using 
PLI, PBI, and our method, USTPL. However, the performance of 
PBI was still lower than PLI and USTPL in these situations. PLI, 
on average, slightly outperformed USTPL when the number of 
trials was large (n = 50, 80), though the difference was not signifi- 
cant. However, USTPL maintained its performance level when the 
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number of trials was decreased to be below the asymptotic limit 
of its significance calculation. This shows that USTPL is robust 
when limited trials are available as opposed to PLI, most likely 
due to the information lost in PLI due to its capturing of only sign 
changes weighted by signal magnitudes whereas USTPL utilizes 
the distribution of phase angles for differentiating the baseline 
from the test population. However, there is a slight advantage 
to using PLI when the number of trials is large and, like in our 
controlled experiment, the baseline contains phase locking due to 
instantaneous effects, such as crosstalk. 

In situations where the prestimulus period is random, while 
the post-stimulus period is phase-locked, we would expect that 
PLV, a parametric method, would have more statistical power 
than USTPL, a non-parametric method, since the assumption of 
a uniform phase distribution holds. However, this is often not the 
case since, due to source modeling or closely spaced sources, any 
correlations would break the assumption of phase distribution 
uniformity. 

Another caveat of our proposed method is that we lose infor- 
mation about differences between phase angle populations. For 
example, if we are comparing two stimulus conditions with 
USTPL, this test alone does not tell us which condition had 
stronger phase locking, or whether one condition leads to a 
greater phase lag. When comparing against a baseline, we assume 
that there is no phase locking within the baseline period. If there 
is suspected phase locking in the baseline, other measures may be 
necessary to supplement differences found in USTPL. However, 



the USTPL method will provide knowledge of those frequency 
and time points amongst ROIs that have significantly different 
phase locking and will thus provide a smaller set of connections 
to investigate with more scrutiny. 

Due to the bivariate nature of phase locking analysis the four 
measures discussed in this paper do not differentiate between 
direct communication between two ROIs and having a third 
source driving both. Methods have been developed to con- 
sider the multivariate network and finding partial phase locking 
(Schelter et al, 2006; Cadieu and Koepsell, 2010; Canolty et al, 
2012). In the future we plan to extend the USTPL model by 
incorporating these multivariate methods. 
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